Spotify Data Analysis and Descriptive Statistics

1. Data Set Description For this project, our team decided to take a data-driven approach to evaluate music. We are using the Spotify Dataset from 1921 to 2020 that consists of over 160,000 different tracks. In order to stay consistent in the evaluation, all data was sourced from the Spotify Web API. The following set of data is combination of Primary, Numerical, Dummy(binary), and Categorical data. This allows us to explore different types of models and draw different insights. Here are all the variables that is included in the data set:

Primary - id: this an unique key comprised of numbers and characters that is assigned to each track generated by Spotify

Numerical - acousticness: range from 0(LOW) to 1(HIGH) - danceability: range from 0(LOW) to 1(HIGH) - energy: range from 0(LOW) to 1(HIGH) - duration_ms: majority range from 200,000 to 300,000 - instrumentalness: range from 0(LOW) to 1(HIGH) - valence: range from 0(LOW) to 1(HIGH) - popularity: range from 0(LOW) to 100(HIGH)* - tempo: majority range from 50(LOW) to 150(high) - liveness: range from 0(LOW) to 1(HIGH) - loudness majority range from -60 to 0 - speechiness: range from 0(LOW) to 1(HIGH) - year: range from 1921 to 2020

Dummy - mode: 0 represents minor and 1 represents major - explicit: 0 represents no explicit content and 1 represents explicit content

Categorical - key: this consists of all different music keys on onctave encoded from 0 to 11 (i.e. C = 0, C# = 1, etc…) - artists: the artist of the track - release_date: the date of release in yyyy-mm-dd format - name: the name of the track

*With our approach on evaluating different tracks, we decided to use popularity as our main dependant variable.

options(repr.matrix.max.rows=100, repr.matrix.max.cols=20)

Import all Libraries

# Libraries
options(warn=-1)
library(ggplot2)
library(dplyr)

Attaching package: ‘dplyr’

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

library(plotly)

Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     
Registered S3 methods overwritten by 'htmltools':
  method               from         
  print.html           tools:rstudio
  print.shiny.tag      tools:rstudio
  print.shiny.tag.list tools:rstudio
Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio

Attaching package: ‘plotly’

The following object is masked from ‘package:ggplot2’:

    last_plot

The following object is masked from ‘package:stats’:

    filter

The following object is masked from ‘package:graphics’:

    layout

library(hrbrthemes)

NOTE: Either Arial Narrow or Roboto Condensed fonts are required to use these themes.
      Please use hrbrthemes::import_roboto_condensed() to install Roboto Condensed and
      if Arial Narrow is not on your system, please see https://bit.ly/arialnarrow

library(forecast)

Registered S3 method overwritten by 'quantmod':
  method            from
  as.zoo.data.frame zoo 

library(xts)

Loading required package: zoo

Attaching package: ‘zoo’

The following objects are masked from ‘package:base’:

    as.Date, as.Date.numeric


Attaching package: ‘xts’

The following objects are masked from ‘package:dplyr’:

    first, last

library(Metrics)

Attaching package: ‘Metrics’

The following object is masked from ‘package:forecast’:

    accuracy

library(psych)

Attaching package: ‘psych’

The following objects are masked from ‘package:ggplot2’:

    %+%, alpha

library(dygraphs)
library(GGally)

Registered S3 method overwritten by 'GGally':
  method from   
  +.gg   ggplot2

library(tidyverse)

Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
── Attaching packages ────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──
✓ tibble  3.0.3     ✓ purrr   0.3.4
✓ tidyr   1.1.2     ✓ stringr 1.4.0
✓ readr   1.4.0     ✓ forcats 0.5.0
── Conflicts ───────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
x psych::%+%()     masks ggplot2::%+%()
x psych::alpha()   masks ggplot2::alpha()
x plotly::filter() masks dplyr::filter(), stats::filter()
x xts::first()     masks dplyr::first()
x dplyr::lag()     masks stats::lag()
x xts::last()      masks dplyr::last()

library(tidyquant)  

Loading required package: lubridate

Attaching package: ‘lubridate’

The following objects are masked from ‘package:base’:

    date, intersect, setdiff, union

Loading required package: PerformanceAnalytics

Attaching package: ‘PerformanceAnalytics’

The following object is masked from ‘package:graphics’:

    legend

Loading required package: quantmod
Loading required package: TTR
Version 0.4-0 included new data defaults. See ?getSymbols.
══ Need to Learn tidyquant? ═══════════════════════════════════════════════════════════════════════════════════════════════════════════════
Business Science offers a 1-hour course - Learning Lab #9: Performance Analysis & Portfolio Optimization with tidyquant!
</> Learn more at: https://university.business-science.io/p/learning-labs-pro </>

library(cranlogs)   
library(corrr)      
library(cowplot)  

Attaching package: ‘cowplot’

The following object is masked from ‘package:lubridate’:

    stamp

library(Metrics)

Initial data visualization

Plotting annual median values for each feature overall. Standardized such that each feature is on the same scale.

df <- read.csv("./data/data.csv")
q1s = data.frame(matrix(ncol=3, nrow=0))
medians = data.frame(matrix(ncol=3, nrow=0))
q3s = data.frame(matrix(ncol=3, nrow=0))
feature_names = sort(c("acousticness", "danceability", "instrumentalness", "energy",
                  "duration_ms", "valence", "tempo", "liveness", "loudness", "speechiness"))
for(i in 1921:2020){
  annual_val = subset(df, year==i)

  for(name in feature_names){
    vals = quantile(annual_val[,name])[2:4]
    q1s <- rbind(q1s, c(i, name, vals[1]))
    medians <- rbind(medians, c(i, name, vals[2]))
    q3s <- rbind(q3s, c(i, name, vals[3]))
  }
}

colnames(q1s) <-c("Year", "Feature", "Value")
colnames(medians) <-c("Year", "Feature", "Value")
colnames(q3s) <-c("Year", "Feature", "Value")

split_q1s <- split(q1s, q1s$Feature)
split_medians <- split(medians, medians$Feature)
split_q3s <- split(q3s, q3s$Feature)

feature_quantiles = array(0, c(100,3,length(feature_names)))

for(i in 1:length(feature_names)){
  for(j in 1:100){
    feature_quantiles[j,1,i] = split_q1s[[i]][[3]][j]
    feature_quantiles[j,2,i] = split_medians[[i]][[3]][j]
    feature_quantiles[j,3,i] = split_q3s[[i]][[3]][j]
  }
}

dimnames(feature_quantiles) <- list(1921:2020,
                                    c("Q1", "Median", "Q3"),
                                    feature_names)



plot(1930,1, xlim = c(1921,2020), ylim = range(0,15), xlab="Year", ylab="Range")
for(i in 1:length(feature_names)){
  name = feature_names[i]
  medians = as.numeric(feature_quantiles[,,name][,"Median"])

  avg = mean(medians)
  medians = medians / avg
  
  lines(x=1921:2020, y=medians, col=i)
}

legend("topright", lty=1, col=c(1,2,3,4,5,6,7,8,9,10), legend=feature_names)

Plotting annual median values for the feature for the top 10% of songs by popularity. Standardized such that each feature is on the same scale.

q1s = data.frame(matrix(ncol=3, nrow=0))
medians = data.frame(matrix(ncol=3, nrow=0))
q3s = data.frame(matrix(ncol=3, nrow=0))
feature_names = sort(c("acousticness", "danceability", "instrumentalness", "energy",
                  "duration_ms", "valence", "tempo", "liveness", "loudness", "speechiness"))
for(i in 1921:2020){
  annual_val = subset(df, year==i)
  
  
  top_10 = floor(nrow(annual_val)*0.1)
  annual_val = annual_val[order(annual_val$popularity, decreasing=TRUE),]
  annual_val = annual_val[1:top_10,]

  for(name in feature_names){
    vals = quantile(annual_val[,name])[2:4]
    q1s <- rbind(q1s, c(i, name, vals[1]))
    medians <- rbind(medians, c(i, name, vals[2]))
    q3s <- rbind(q3s, c(i, name, vals[3]))
  }
}

colnames(q1s) <-c("Year", "Feature", "Value")
colnames(medians) <-c("Year", "Feature", "Value")
colnames(q3s) <-c("Year", "Feature", "Value")

split_q1s <- split(q1s, q1s$Feature)
split_medians <- split(medians, medians$Feature)
split_q3s <- split(q3s, q3s$Feature)

feature_quantiles = array(0, c(100,3,length(feature_names)))

for(i in 1:length(feature_names)){
  for(j in 1:100){
    feature_quantiles[j,1,i] = split_q1s[[i]][[3]][j]
    feature_quantiles[j,2,i] = split_medians[[i]][[3]][j]
    feature_quantiles[j,3,i] = split_q3s[[i]][[3]][j]
  }
}

dimnames(feature_quantiles) <- list(1921:2020,
                                    c("Q1", "Median", "Q3"),
                                    feature_names)


plot(1930,1, xlim = c(1921,2020), ylim = range(0,17), xlab="Year", ylab="Range")
for(i in 1:length(feature_names)){
  name = feature_names[i]
  medians = as.numeric(feature_quantiles[,,name][,"Median"])
  avg = mean(medians)
  medians = medians / avg
  
  lines(x=1921:2020, y=medians, col=i)
}

legend("topright", lty=1, col=c(1,2,3,4,5,6,7,8,9,10), legend=feature_names)

####Insights:
Both these graphs appear to have similar shapes, with instrumentalness being the most widely varying, and the rest varying marginally. However, for the speechiness variable,

Now plotting absolute values for each feature as the median values for overall, as well as just the most popular songs.

MAPE is the difference between the most popular median and the overall median; R is the correlation between overall median and most popular median.

These values will determine whether or not the overall median is a useful predictor of what the median for popular songs will be. This is important because it will be easier to forecast overall values into the future, given they are much smoother, and these forecast overall values can be used to predict the most popular feature value.

# Get the list of features I need to work with

library(Metrics)

q1s = data.frame(matrix(ncol=3, nrow=0))
medians = data.frame(matrix(ncol=3, nrow=0))
best_medians = data.frame(matrix(ncol=3, nrow=0))
q3s = data.frame(matrix(ncol=3, nrow=0))
feature_names = sort(c("acousticness", "danceability", "instrumentalness", "energy",
                  "duration_ms", "valence", "tempo", "liveness", "loudness", "speechiness"))
for(i in 1921:2020){
  annual_val = subset(df, year==i)
  top_10 = floor(nrow(annual_val)*0.1)
  annual_val = annual_val[order(annual_val$popularity, decreasing=TRUE),]
  best_ones = annual_val[1:top_10,]

  for(name in feature_names){
    vals = quantile(annual_val[,name])[2:4]
    best_vals = quantile(best_ones[,name])[2:4]
    best_medians <- rbind(best_medians, c(i, name, best_vals[2]))
    medians <- rbind(medians, c(i, name, vals[2]))
  }
}

colnames(medians) <-c("Year", "Feature", "Value")
colnames(best_medians) <- c("Year", "Feature", "Value")

split_medians <- split(medians, medians$Feature)
split_best_medians <- split(best_medians, best_medians$Feature)

feature_quantiles = matrix(0, ncol=length(feature_names), nrow=100)
best_feature_quantiles = matrix(0, ncol=length(feature_names), nrow=100)

for(i in 1:length(feature_names)){
  for(j in 1:100){
    feature_quantiles[j,i] <- split_medians[[i]][[3]][j]
    best_feature_quantiles[j,i] = split_best_medians[[i]][[3]][j]
  }
}

dimnames(feature_quantiles) <- list(1921:2020,
                                    feature_names)
dimnames(best_feature_quantiles) <- list(1921:2020,
                                    feature_names)




for(i in 1:length(feature_names)){
  name = feature_names[i]
  scores = as.numeric(feature_quantiles[,i])
  best_scores = as.numeric(best_feature_quantiles[,i])
  mape = mape(best_scores, scores)
  mape = round(mape * 100,3)
  r  = cor(scores, best_scores)
  pval = t.test(scores, )
  # Normalize it????? / make it a function of its own mean
  
  smallest = min(min(scores), min(best_scores))
  biggest = max(max(scores), max(best_scores))
  
  plot(x=1921:2020, y=scores, col=1, type="l", ylim=c(smallest, biggest), main=name, sub=paste("MAPE: ", mape, "%;", "R:", r), xlab="Year", ylab="Median Value")
  lines(x=1921:2020, y=best_scores, col=2)
  legend("topright", lty=1, col=c(1,2), legend=c("Overall", "Most Popular"))
  title()
}

Insights:

Acousticness: A high MAPE when comparing the most popular songs with the overall trends indicates that there is a large difference between the two (likely caused in the 1930s-1950s). The R value is high which means there is a high correlation between the most popular songs and the overall songs. Because acousticness seems to impact popular songs and there is a large change over time, we will move forward with predicting the acousticness of songs in the future. Danceability: A low MAPE indicates not a lot of difference between popular songs and overall songs danceability. However, a high R value means there is a high correlation between the most popular songs and the overall songs. Because danceability seems to impact popular songs (due to the high correlation) and there is a some change over time, we will move forward with predicting the danceability of songs in the future. Duration (ms): There is a low MAPE which means that popular songs have about the same length as all other songs. The R value is relatively high which means there is a high correlation between the most popular songs and the overall songs. Because duration seems to impact popular songs (due to the high correlation) and there is a some change over time, we will move forward with predicting the duration of songs in the future. Energy: A high MAPE when comparing the most popular songs with the overall trends indicates that there is a large difference between the two. The R value is high which means there is a high correlation between the most popular songs and the overall songs. Because energy seems to impact popular songs (due to the high MAPE and correlation) and there are large changes over time, we will move forward with predicting the energy of songs in the future. Instrumentalness: The popular songs and overall music trends are pretty much the exact same after the 1950s, indicating that instrumentalness does not affect popularity. Because instrumentalness does not impact popular songs and there is a flatline after the 1950s (indicating no change over time), we will not move forward with predicting the instrumentalness of songs in the future. Liveness: A high MAPE when comparing the most popular songs with the overall trends indicates that there is a large difference between the two. The R value is high which means there is a high correlation between the most popular songs and the overall songs. Because liveness seems to impact popular songs (due to the high MAPE and correlation) and there is a some change over time, we will move forward with predicting the liveness of songs in the future. Loudness: A high MAPE when comparing the most popular songs with the overall trends indicates that there is a large difference between the two. The R value is high which means there is a high correlation between the most popular songs and the overall songs. Because loudness seems to impact popular songs (due to the high MAPE and correlation) and there is a significant change over time, we will move forward with predicting the loudness of songs in the future. Speechiness: A high MAPE when comparing the most popular songs with the overall trends indicates that there is a large difference between the two. It’s important to note that the driver of the high MAPE is due to the two large spikes between 1930 and 1950, while speechiness levels between popular songs and the overall industry in the last 70 years have generally been the same. The R value is low which means there is a low correlation between the most popular songs and the overall songs. We will not move forward with predicting speechiness, due to the lack of change in speechiness over time and low R. Tempo: The MAPE is low, inidcating that there isn’t a significant difference in tempo between the popular songs and the overall music trends. The R is relatively high, so there is possibly a correlation between popular songs and the overall songs. Because tempo seems to impact popular songs (due to the high correlation) and there is a some change over time, we will move forward with predicting the tempo of songs in the future. Valence: A high MAPE when comparing the most popular songs with the overall trends indicates that there is a large difference between the two in terms of valence. The R value is high which means there is a high correlation between the most popular songs and the overall songs. Because valence seems to impact popular songs (due to the high MAPE and high correlation) and there is a some change over time, we will move forward with predicting the valence of songs in the future.

Train time series models to forecast the future models.

# Set up the libraries and the training/testing amounts
library("smooth")

Loading required package: greybox
Package "greybox", v0.6.3 loaded.


Attaching package: ‘greybox’

The following object is masked from ‘package:lubridate’:

    hm

The following object is masked from ‘package:tidyr’:

    spread

This is package "smooth", v2.6.0

Attaching package: ‘smooth’

The following object is masked from ‘package:TTR’:

    lags

library("forecast")
library("nnfor")
training.percent = 0.95
nTrain = 100*training.percent
nTest = 100*(1-training.percent)

Accousticness Models

acousticness = ts(as.numeric(feature_quantiles[,"acousticness"]) , start=1921)
acousticness.train = subset(acousticness, start=1, end=nTrain)
acousticness.test = subset(acousticness, start = (nTrain+1), end =(nTrain+nTest))

#ses
acousticness.train.ses <- ses(acousticness.train, h=nTest)
acousticness.ses.mape = mape(acousticness.train.ses$mean, acousticness.test)
plot(acousticness.train.ses, main=paste("Acousticness", "SES"), sub=paste("MAPE:", round(acousticness.ses.mape*100, 3), "%"), xlab="Year", ylab="Median Value")
lines(acousticness)

#SARIMA
acousticness.sarima.model <- auto.arima(acousticness.train)
acousticness.sarima <- forecast(acousticness.sarima.model, h=nTest)
acousticness.sarima.mape = mape(acousticness.sarima$mean, acousticness.test)
plot(acousticness.sarima, main=paste("Acousticness", "SARIMA"), sub=paste("MAPE:", round(acousticness.sarima.mape*100,3), "%"), xlab="Year", ylab="Median Value")
lines(acousticness)

#Neural Network
acousticness.train.nn <- elm(acousticness.train)
acousticness.nn.forecast <- forecast(acousticness.train.nn, h=nTest)
acousticness.nn.mape = mape(acousticness.nn.forecast$mean, acousticness.test)
plot(acousticness, main=paste("Acousticness", "NN"), sub=paste("MAPE:", round(acousticness.nn.mape*100, 3), "%"), xlab="Year", ylab="Median Value")
lines(ts(acousticness.nn.forecast$mean, start=1921+nTrain), col=2)

Here we see that the best method to predict future acousticness median values for the overall data set is through SES. However, given that this best MAPE is roughly 37%, this feature cannot be reliably forecast into the future.

Energy models

energy = ts(as.numeric(feature_quantiles[,"energy"]) , start=1921)
energy.train = subset(energy, start=1, end=nTrain)
energy.test = subset(energy, start = (nTrain+1), end =(nTrain+nTest))

#SES
energy.train.ses <- ses(energy.train, h=nTest)
energy.ses.mape = mape(energy.train.ses$mean, energy.test)
plot(energy.train.ses, main=paste("Energy", "SES"), sub=paste("MAPE:", round(energy.ses.mape* 100, 3), "%"), xlab="Year", ylab="Median Value")
lines(energy)

#SARIMA
energy.sarima.model <- auto.arima(energy.train)
energy.sarima <- forecast(energy.sarima.model, h=nTest)
energy.sarima.mape = mape(energy.sarima$mean, energy.test)
plot(energy.sarima, main=paste("Energy", "SARIMA"), sub=paste("MAPE:", round(energy.sarima.mape* 100, 3), "%"), xlab="Year", ylab="Median Value")
lines(energy)

#Neural Network
energy.train.nn <- elm(energy.train)
energy.nn.forecast <- forecast(energy.train.nn, h=nTest)
energy.nn.mape = mape(energy.nn.forecast$mean, energy.test)
plot(energy, main=paste("Energy", "NN"), sub=paste("MAPE:", round(energy.nn.mape* 100, 3), "%"), xlab="Year", ylab="Median Value")
lines(ts(energy.nn.forecast$mean, start=1921+nTrain), col=2)

Here we see that the best method to predict future energy median values for the overall data set is through SES. Given that this best MAPE is roughly 8%, this feature can be fairly reliably forecast into the future.

Danceability models

danceability = ts(as.numeric(feature_quantiles[,"danceability"]) , start=1921)
danceability.train = subset(danceability, start=1, end=nTrain)
danceability.test = subset(danceability, start = (nTrain+1), end =(nTrain+nTest))

#SES
danceability.train.ses <- ses(danceability.train, h=nTest)
danceability.ses.mape = mape(danceability.train.ses$mean, danceability.test)
plot(danceability.train.ses, main=paste("danceability", "SES"), sub=paste("MAPE:", round(danceability.ses.mape* 100, 3), "%"), xlab="Year", ylab="Median Value")
lines(danceability)

#SARIMA
danceability.sarima.model <- auto.arima(danceability.train)
danceability.sarima <- forecast(danceability.sarima.model, h=nTest)
danceability.sarima.mape = mape(danceability.sarima$mean, danceability.test)
plot(danceability.sarima, main=paste("danceability", "SARIMA"), sub=paste("MAPE:", round(danceability.sarima.mape* 100, 3), "%"), xlab="Year", ylab="Median Value")
lines(danceability)

#Neural Network
danceability.train.nn <- elm(danceability.train)
danceability.nn.forecast <- forecast(danceability.train.nn, h=nTest)
danceability.nn.mape = mape(danceability.nn.forecast$mean, danceability.test)
plot(danceability, main=paste("danceability", "NN"), sub=paste("MAPE:", round(danceability.nn.mape* 100, 3), "%"), xlab="Year", ylab="Median Value")
lines(ts(danceability.nn.forecast$mean, start=1921+nTrain), col=2)

Here we see that the best method to predict future danceability median values for the overall data set is through NN. However, given that this best MAPE is roughly 10.4%, this feature cannot be reliably forecast into the future.

Duration models

duration = ts(as.numeric(feature_quantiles[,"duration_ms"]) , start=1921)
duration.train = subset(duration, start=1, end=nTrain)
duration.test = subset(duration, start = (nTrain+1), end =(nTrain+nTest))

#SES
duration.train.ses <- ses(duration.train, h=nTest)
duration.ses.mape = mape(duration.train.ses$mean, duration.test)
plot(duration.train.ses, main=paste("Duration", "SES"), sub=paste("MAPE:",round(duration.ses.mape*100, 3), "%"), xlab="Year", ylab="Median Value")
lines(duration)

#SARIMA
duration.sarima.model <- auto.arima(duration.train)
duration.sarima <- forecast(duration.sarima.model, h=nTest)
duration.sarima.mape = mape(duration.sarima$mean, duration.test)
plot(duration.sarima, main=paste("Duration", "SARIMA"), sub=paste("MAPE:", round(duration.sarima.mape*100, 3), "%"), xlab="Year", ylab="Median Value")
lines(duration)

#Neural Network
duration.train.nn <- elm(duration.train)
duration.nn.forecast <- forecast(duration.train.nn, h=nTest)
duration.nn.mape = mape(duration.nn.forecast$mean, duration.test)
plot(duration, main=paste("Duration", "NN"), sub=paste("MAPE:", round(duration.nn.mape*100, 3), "%"), xlab="Year", ylab="Median Value")
lines(ts(duration.nn.forecast$mean, start=1921+nTrain), col=2)

Here we see that the best method to predict future duration median values for the overall data set is through NN. Given that this best MAPE is roughly 8%, this feature can be fairly reliably forecast into the future.

Valence

valence = ts(as.numeric(feature_quantiles[,"valence"]) , start=1921)
valence.train = subset(valence, start=1, end=nTrain)
valence.test = subset(valence, start = (nTrain+1), end =(nTrain+nTest))

#SES
valence.train.ses <- ses(valence.train, h=nTest)
valence.ses.mape = mape(valence.train.ses$mean, valence.test)
plot(valence.train.ses, main=paste("valence", "SES"), sub=paste("MAPE:",round(valence.ses.mape*100, 3), "%"), xlab="Year", ylab="Median Value")
lines(valence)

#SARIMA
valence.sarima.model <- auto.arima(valence.train)
valence.sarima <- forecast(valence.sarima.model, h=nTest)
valence.sarima.mape = mape(valence.sarima$mean, valence.test)
plot(valence.sarima, main=paste("valence", "SARIMA"), sub=paste("MAPE:", round(valence.sarima.mape*100, 3), "%"), xlab="Year", ylab="Median Value")
lines(valence)

#Neural Network
valence.train.nn <- elm(valence.train)
valence.nn.forecast <- forecast(valence.train.nn, h=nTest)
valence.nn.mape = mape(valence.nn.forecast$mean, valence.test)
plot(valence, main=paste("valence", "NN"), sub=paste("MAPE:", round(valence.nn.mape*100, 3), "%"), xlab="Year", ylab="Median Value")
lines(ts(valence.nn.forecast$mean, start=1921+nTrain), col=2)

Insights:

Relatively lower MAPE indicates that we can predict valence over time. The best model to use is the SARIMA model because it has the lowest MAPE at 6.7%.

Tempo model

tempo = ts(as.numeric(feature_quantiles[,"tempo"]) , start=1921)
tempo.train = subset(tempo, start=1, end=nTrain)
tempo.test = subset(tempo, start = (nTrain+1), end =(nTrain+nTest))

#SES
tempo.train.ses <- ses(tempo.train, h=nTest)
tempo.ses.mape = mape(tempo.train.ses$mean, tempo.test)
plot(tempo.train.ses, main=paste("tempo", "SES"), sub=paste("MAPE:",round(tempo.ses.mape*100, 3), "%"), xlab="Year", ylab="Median Value")
lines(tempo)

#SARIMA
tempo.sarima.model <- auto.arima(tempo.train)
tempo.sarima <- forecast(tempo.sarima.model, h=nTest)
tempo.sarima.mape = mape(tempo.sarima$mean, tempo.test)
plot(tempo.sarima, main=paste("tempo", "SARIMA"), sub=paste("MAPE:", round(tempo.sarima.mape*100, 3), "%"), xlab="Year", ylab="Median Value")
lines(tempo)

#Neural Network
tempo.train.nn <- elm(tempo.train)
tempo.nn.forecast <- forecast(tempo.train.nn, h=nTest)
tempo.nn.mape = mape(tempo.nn.forecast$mean, tempo.test)
plot(tempo, main=paste("tempo", "NN"), sub=paste("MAPE:", round(tempo.nn.mape*100, 3), "%"), xlab="Year", ylab="Median Value")
lines(ts(tempo.nn.forecast$mean, start=1921+nTrain), col=2)

Insights:

Here we see that the best method to predict future tempo median values for the overall data set is through SES. Given that this best MAPE is roughly 1%, this feature can be fairly reliably forecast into the future.

Liveness

liveness = ts(as.numeric(feature_quantiles[,"liveness"]) , start=1921)
liveness.train = subset(liveness, start=1, end=nTrain)
liveness.test = subset(liveness, start = (nTrain+1), end =(nTrain+nTest))

#SES
liveness.train.ses <- ses(liveness.train, h=nTest)
liveness.ses.mape = mape(liveness.train.ses$mean, liveness.test)
plot(liveness.train.ses, main=paste("liveness", "SES"), sub=paste("MAPE:",round(liveness.ses.mape*100, 3), "%"), xlab="Year", ylab="Median Value")
lines(liveness)

#SARIMA
liveness.sarima.model <- auto.arima(liveness.train)
liveness.sarima <- forecast(liveness.sarima.model, h=nTest)
liveness.sarima.mape = mape(liveness.sarima$mean, liveness.test)
plot(liveness.sarima, main=paste("liveness", "SARIMA"), sub=paste("MAPE:", round(liveness.sarima.mape*100, 3), "%"), xlab="Year", ylab="Median Value")
lines(liveness)

#Neural Network
liveness.train.nn <- elm(liveness.train)
liveness.nn.forecast <- forecast(liveness.train.nn, h=nTest)
liveness.nn.mape = mape(liveness.nn.forecast$mean, liveness.test)
plot(liveness, main=paste("liveness", "NN"), sub=paste("MAPE:", round(liveness.nn.mape*100, 3), "%"), xlab="Year", ylab="Median Value")
lines(ts(liveness.nn.forecast$mean, start=1921+nTrain), col=2)

Insights:

Here we see that the best method to predict future liveness median values for the overall data set is through SARIMA Given that this best MAPE is roughly 2.6%, this feature can be fairly reliably forecast into the future.

Loudness

loudness = ts(as.numeric(feature_quantiles[,"loudness"]) , start=1921)
loudness.train = subset(loudness, start=1, end=nTrain)
loudness.test = subset(loudness, start = (nTrain+1), end =(nTrain+nTest))

#SES
loudness.train.ses <- ses(loudness.train, h=nTest)
loudness.ses.mape = mape(loudness.train.ses$mean, loudness.test)
plot(loudness.train.ses, main=paste("loudness", "SES"), sub=paste("MAPE:",round(loudness.ses.mape*100, 3), "%"), xlab="Year", ylab="Median Value")
lines(loudness)

#SARIMA
loudness.sarima.model <- auto.arima(loudness.train)
loudness.sarima <- forecast(loudness.sarima.model, h=nTest)
loudness.sarima.mape = mape(loudness.sarima$mean, loudness.test)
plot(loudness.sarima, main=paste("loudness", "SARIMA"), sub=paste("MAPE:", round(loudness.sarima.mape*100, 3), "%"), xlab="Year", ylab="Median Value")
lines(loudness)

#Neural Network
loudness.train.nn <- elm(loudness.train)
loudness.nn.forecast <- forecast(loudness.train.nn, h=nTest)
loudness.nn.mape = mape(loudness.nn.forecast$mean, loudness.test)
plot(loudness, main=paste("loudness", "NN"), sub=paste("MAPE:", round(loudness.nn.mape*100, 3), "%"), xlab="Year", ylab="Median Value")
lines(ts(loudness.nn.forecast$mean, start=1921+nTrain), col=2)

Insights:

Here we see that the best method to predict future loudness median values for the overall data set is through SES Given that this best MAPE is roughly 5%, this feature can be fairly reliably forecast into the future.

Linear Model for Overall Feature Median–>Most Popular Songs’ Feature Median

We will only consider models that had a low enough MAPEs (10% chosen as a threshold) in their best timeseries forecasting model.

Here MAPE is the difference between the actual median for most popular songs and the predicted median.

#Energy
best_energy = ts(as.numeric(best_feature_quantiles[,"energy"]), start=1921, end=2020) 
best_energy.train = subset(best_energy, start = 1, end = nTrain)
best_energy.test = subset(best_energy, start = nTrain+1, end = nTrain+nTest)
train_df = data.frame(best = best_energy.train, overall = energy.train)

best = best_energy.train
overall = energy.train

energy.lm <- lm(best ~ overall)
test_df <- data.frame(overall = energy.test)
best_energy.test.predict <- predict(energy.lm, newdata=test_df, interval="prediction")
best_energy_test_predict.mape <- mape(best_energy.test[1:length(best_energy.test)], best_energy.test.predict[,1])
plot(best_energy,xlab="Year", ylab="Energy", main=paste("Energy Popularity Prediction"), sub = paste("MAPE: ", round(best_energy_test_predict.mape*100, 3), "%"))
lines(ts(best_energy.test.predict[,1], start = 1921+nTrain+1, end = 2020), col = 2)

#Duration
best_duration = ts(as.numeric(best_feature_quantiles[,"duration_ms"]), start=1921, end=2020) 
best_duration.train = subset(best_duration, start = 1, end = nTrain)
best_duration.test = subset(best_duration, start = nTrain+1, end = nTrain+nTest)
train_df = data.frame(best = best_duration.train, overall = duration.train)

best = best_duration.train
overall = duration.train

duration.lm <- lm(best ~ overall)
test_df <- data.frame(overall = duration.test)
best_duration.test.predict <- predict(duration.lm, newdata=test_df, interval="prediction")
best_duration_test.predict.mape <- mape(best_duration.test[1:length(best_duration.test)], best_duration.test.predict[,1])
plot(best_duration,xlab="Year", ylab="duration", main=paste("duration Popularity Prediction"), sub = paste("MAPE: ", round(best_duration_test.predict.mape*100, 3), "%"))
lines(ts(best_duration.test.predict[,1], start = 1921+nTrain+1, end = 2020), col = 2)

#Accousticness
best_acousticness = ts(as.numeric(best_feature_quantiles[,"acousticness"]), start=1921, end=2020) 
best_acousticness.train = subset(best_acousticness, start = 1, end = nTrain)
best_acousticness.test = subset(best_acousticness, start = nTrain+1, end = nTrain+nTest)
train_df = data.frame(best = best_acousticness.train, overall = acousticness.train)

best = best_acousticness.train
overall = acousticness.train

acousticness.lm <- lm(best ~ overall)
test_df <- data.frame(overall = acousticness.test)
best_acousticness.test.predict <- predict(acousticness.lm, newdata=test_df, interval="prediction")
best_acousticness_test.predict.mape <- mape(best_acousticness.test[1:length(best_acousticness.test)], best_acousticness.test.predict[,1])
plot(best_acousticness,xlab="Year", ylab="acousticness", main=paste("acousticness Popularity Prediction"), sub = paste("MAPE: ", round(best_acousticness_test.predict.mape*100, 3), "%"))
lines(ts(best_acousticness.test.predict[,1], start = 1921+nTrain+1, end = 2020), col = 2)

#Valence
best_valence = ts(as.numeric(best_feature_quantiles[,"valence"]), start=1921, end=2020) 
best_valence.train = subset(best_valence, start = 1, end = nTrain)
best_valence.test = subset(best_valence, start = nTrain+1, end = nTrain+nTest)
train_df = data.frame(best = best_valence.train, overall = valence.train)

best = best_valence.train
overall = valence.train

valence.lm <- lm(best ~ overall)
test_df <- data.frame(overall = valence.test)
best_valence.test.predict <- predict(valence.lm, newdata=test_df, interval="prediction")
best_valence_test.predict.mape <- mape(best_valence.test[1:length(best_valence.test)], best_valence.test.predict[,1])
plot(best_valence,xlab="Year", ylab="Valence", main=paste("Valence Popularity Prediction"), sub = paste("MAPE: ", round(best_valence_test.predict.mape*100, 3), "%"))
lines(ts(best_valence.test.predict[,1], start = 1921+nTrain+1, end = 2020), col = 2)

#Tempo
best_tempo = ts(as.numeric(best_feature_quantiles[,"tempo"]), start=1921, end=2020) 
best_tempo.train = subset(best_tempo, start = 1, end = nTrain)
best_tempo.test = subset(best_tempo, start = nTrain+1, end = nTrain+nTest)
train_df = data.frame(best = best_tempo.train, overall = tempo.train)

best = best_tempo.train
overall = tempo.train

tempo.lm <- lm(best ~ overall)
test_df <- data.frame(overall = tempo.test)
best_tempo.test.predict <- predict(tempo.lm, newdata=test_df, interval="prediction")
best_tempo_test.predict.mape <- mape(best_tempo.test[1:length(best_tempo.test)], best_tempo.test.predict[,1])
plot(best_tempo,xlab="Year", ylab="tempo", main=paste("tempo Popularity Prediction"), sub = paste("MAPE: ", round(best_tempo_test.predict.mape*100, 3), "%"))
lines(ts(best_tempo.test.predict[,1], start = 1921+nTrain+1, end = 2020), col = 2)

#Liveness
best_liveness = ts(as.numeric(best_feature_quantiles[,"liveness"]), start=1921, end=2020) 
best_liveness.train = subset(best_liveness, start = 1, end = nTrain)
best_liveness.test = subset(best_liveness, start = nTrain+1, end = nTrain+nTest)
train_df = data.frame(best = best_liveness.train, overall = liveness.train)

best = best_liveness.train
overall = liveness.train

liveness.lm <- lm(best ~ overall)
test_df <- data.frame(overall = liveness.test)
best_liveness.test.predict <- predict(liveness.lm, newdata=test_df, interval="prediction")
best_liveness_test.predict.mape <- mape(best_liveness.test[1:length(best_liveness.test)], best_liveness.test.predict[,1])
plot(best_liveness,xlab="Year", ylab="liveness", main=paste("liveness Popularity Prediction"), sub = paste("MAPE: ", round(best_liveness_test.predict.mape*100, 3), "%"))
lines(ts(best_liveness.test.predict[,1], start = 1921+nTrain+1, end = 2020), col = 2)

#Loudness
best_loudness = ts(as.numeric(best_feature_quantiles[,"loudness"]), start=1921, end=2020) 
best_loudness.train = subset(best_loudness, start = 1, end = nTrain)
best_loudness.test = subset(best_loudness, start = nTrain+1, end = nTrain+nTest)
train_df = data.frame(best = best_loudness.train, overall = loudness.train)

best = best_loudness.train
overall = loudness.train

loudness.lm <- lm(best ~ overall)
test_df <- data.frame(overall = loudness.test)
best_loudness.test.predict <- predict(loudness.lm, newdata=test_df, interval="prediction")
best_loudness_test.predict.mape <- mape(best_loudness.test[1:length(best_loudness.test)], best_loudness.test.predict[,1])
plot(best_loudness,xlab="Year", ylab="loudness", main=paste("loudness Popularity Prediction"), sub = paste("MAPE: ", round(best_loudness_test.predict.mape*100, 3), "%"))
lines(ts(best_loudness.test.predict[,1], start = 1921+nTrain+1, end = 2020), col = 2)

Insights:

Energy: A low MAPE around 4% indicates that we can use the overall trend to identify the energy of popular songs Duration: A low MAPE around 1% indicates that we can use the overall trend to identify the duration of popular songs Acousticness: A high MAPE around 16% indicates that we can not use the overall trend to identify the acousticness of popular songs. We will not include predictions for acousticness of popular songs. Valence: A low MAPE around 4.6% indicates that we can use the overall trend to identify the valence of popular songs Tempo: A low MAPE around 2% indicates that we can use the overall trend to identify the tempo of popular songs Liveness: A low MAPE around 4.3% indicates that we can use the overall trend to identify the liveness of popular songs Loudness: A low MAPE around 5% indicates that we can use the overall trend to identify the loudness of popular songs

Predict features for the overall market trends for the next five years.

Only including features whose trained linear models had a MAPE below 10%, as these are the features whose most popular medians will be most reliably predicted by their forecast overall values.

Using the best timeseries method seen for each feature.

Energy: SES Model

energy.ses <- ses(energy, h=nTest)
plot(energy.ses, main=paste("Energy", "SES"), xlab="Year", ylab="Median Value")
lines(energy)

Duration: NN Model

duration.nn <- elm(duration)
duration.nn.predict <- forecast(duration.nn, h=nTest)
plot(duration.nn.predict, main=paste("Duration", "NN"), xlab="Year", ylab="Median Value")

Tempo: SES Model

tempo.ses <- ses(tempo, h=nTest)
plot(tempo.ses, main=paste("Tempo", "SES"), xlab="Year", ylab="Median Value")
lines(tempo)

Valence: SARIMA Model

valence.sarima <- auto.arima(valence)
valence.sarima.predict <- forecast(valence.sarima, h=nTest)
plot(valence.sarima.predict, main=paste("Valence", "SARIMA"), xlab="Year", ylab="Median Value")
lines(valence)

Liveness: SARIMA Model

liveness.sarima <- auto.arima(liveness)
liveness.sarima.predict <- forecast(liveness.sarima, h=nTest)
plot(liveness.sarima.predict, main=paste("liveness", "SARIMA"), xlab="Year", ylab="Median Value")
lines(liveness)

Loudness: SES Model

loudness.ses <- ses(loudness, h=nTest)
plot(loudness.ses, main=paste("loudness", "SES"), xlab="Year", ylab="Median Value")
lines(loudness)

Now predicting the median values for the most popular songs for this set of features.

The forecast overall median values and linear model predicting the most popular median values for each feature given overall median value will be used to predict the future median values for popular songs for each feature.

Energy

best = best_energy
overall = energy


energy.lm <- lm(best ~ overall)
predict_df <- data.frame(overall = energy.ses$mean)
energy.predict <- predict(energy.lm, newdata=predict_df, interval="prediction")

smallest = min(min(best), min(overall))
biggest = max(max(best), max(overall))
plot(energy, xlab="Year", ylab="Median Value", ylim=c(smallest, biggest), main="Most Popular Song's Energy Feature Prediction", sub=paste("Projected median feature value for most popular songs:", paste(round(energy.ses$mean, 3), collapse=",") ))
lines(best_energy, col = 2)

lines(ts(energy.predict[,1], start=2020+1, end=2020+nTest), col = 3)
legend("bottomright", lty=1, col=c(1,2,3), legend=c("Overall", "Most Popular", "Most Popular Projection"))

Duration

best = best_duration
overall = duration

duration.lm <- lm(best ~ overall)
predict_df <- data.frame(overall = duration.nn.predict$mean)
duration.predict <- predict(duration.lm, newdata=predict_df, interval="prediction")

smallest = min(min(best), min(overall))
biggest = max(max(best), max(overall))
plot(duration, xlab="Year", ylab="Median Value", ylim=c(smallest, biggest), main="Most Popular Song's Duration Feature Prediction", sub=paste("Projected median feature value for most popular songs:", paste(round(duration.nn.predict$mean), collapse=",") ))
lines(best_duration, col = 2)

lines(ts(duration.predict[,1], start=2020+1, end=2020+nTest), col = 3)
legend("topright", lty=1, col=c(1,2,3), legend=c("Overall", "Most Popular", "Most Popular Projection"))

Tempo

best = best_tempo
overall = tempo

tempo.lm <- lm(best ~ overall)
predict_df <- data.frame(overall = tempo.ses$mean)
tempo.predict <- predict(tempo.lm, newdata=predict_df, interval="prediction")


smallest = min(min(best), min(overall))
biggest = max(max(best), max(overall))
plot(tempo, xlab="Year", ylab="Median Value", ylim=c(smallest, biggest), main="Most Popular Song's Tempo Feature Prediction", sub=paste("Projected median feature value for most popular songs:", paste(round(tempo.ses$mean), collapse=",") ))
lines(best_tempo, col = 2)

lines(ts(tempo.predict[,1], start=2020+1, end=2020+nTest), col = 3)
legend("bottomright", lty=1, col=c(1,2,3), legend=c("Overall", "Most Popular", "Most Popular Projection"))

Valence

best = best_valence
overall = valence

valence.lm <- lm(best ~ overall)
predict_df <- data.frame(overall = valence.sarima.predict$mean)
valence.predict <- predict(valence.lm, newdata=predict_df, interval="prediction")


smallest = min(min(best), min(overall))
biggest = max(max(best), max(overall))
plot(valence, xlab="Year", ylab="Median Value", ylim=c(smallest, biggest), main="Most Popular Song's valence Feature Prediction", sub=paste("Projected median feature value for most popular songs:", paste(round(valence.sarima.predict$mean,3), collapse=",") ))
lines(best_valence, col = 2)

lines(ts(valence.predict[,1], start=2020+1, end=2020+nTest), col = 3)
legend("bottomright", lty=1, col=c(1,2,3), legend=c("Overall", "Most Popular", "Most Popular Projection"))

Liveness

best = best_liveness
overall = liveness

liveness.lm <- lm(best ~ overall)
predict_df <- data.frame(overall = liveness.sarima.predict$mean)
liveness.predict <- predict(liveness.lm, newdata=predict_df, interval="prediction")

smallest = min(min(best), min(overall))
biggest = max(max(best), max(overall))
plot(liveness, xlab="Year", ylab="Median Value", ylim=c(smallest, biggest), main="Most Popular Song's Liveness Feature Prediction", sub=paste("Projected median feature value for most popular songs:", paste(round(liveness.sarima.predict$mean,3), collapse=",") ))
lines(best_liveness, col = 2)

lines(ts(liveness.predict[,1], start=2020+1, end=2020+nTest), col = 3)
legend("topright", lty=1, col=c(1,2,3), legend=c("Overall", "Most Popular", "Most Popular Projection"))

Loudness

best = best_loudness
overall = loudness

loudness.lm <- lm(best ~ overall)
predict_df <- data.frame(overall = loudness.ses$mean)
loudness.predict <- predict(loudness.lm, newdata=predict_df, interval="prediction")

smallest = min(min(best), min(overall))
biggest = max(max(best), max(overall))
plot(loudness, xlab="Year", ylab="Median Value", ylim=c(smallest, biggest), main="Most Popular Song's Loudness Feature Prediction", sub=paste("Projected median feature value for most popular songs:", paste(round(loudness.ses$mean,1), collapse=",") ))
lines(best_loudness, col = 2)

lines(ts(loudness.predict[,1], start=2020+1, end=2020+nTest), col = 3)
legend("bottomright", lty=1, col=c(1,2,3), legend=c("Overall", "Most Popular", "Most Popular Projection"))

Implications:

Using these values above for the predicted median value for the most popular songs, for each feature, music curators and record labels can predict which “sounds” will be most popular in the future. This information can then be used for future decision making with respect to which artists to sign and which sounds to push.